'''
Nimplanted = 28.
NscatterDisk = 1e4
RatioSDtotal = 0.0529
RatioSDtotal_unc = 0.0034
PimplSD = Nimplanted/NscatterDisk * RatioSDtotal

PearthImpl = 0.12
PearthImpl_unc = 0.04
MassBelt = 4e-4
MassLV = 4.8e-3
MassLV_unc = 1.63e-3
Nasteroids = 607

Nplanetsimals = Nasteroids * MassLV / MassBelt / PearthImpl
Nimpl = PimplSD * Nplanetsimals

chi=Nimplanted/NscatterDisk
Nimpl_unc = chi *  Nasteroids / MassBelt * MassLV / PearthImpl * RatioSDtotal_unc + \
            chi *  RatioSDtotal * Nasteroids / MassBelt * 1 / PearthImpl * MassLV_unc + \
            chi *  RatioSDtotal * Nasteroids / MassBelt * MassLV / PearthImpl**2 * PearthImpl_unc

print(Nimpl, Nimpl_unc)
'''
import argparse
parser = argparse.ArgumentParser(description='Process and save images from a camera')
parser.add_argument("-f","--full",action="store_true",help="entire IMB")
args = parser.parse_args()

#### INNER MAIN BELT
if (args.full==True):
    nImpSim = 28.
    nImpSim_unc = nImpSim**0.5
    nSdSim = 1.e4
    sd_to_total = 0.053
    sd_to_total_unc = 0.003
else:
    nImpSim = 6.
    nImpSim_unc = nImpSim**0.5
    nSdSim = 1.e4
    sd_to_total = 0.053
    sd_to_total_unc = 0.003

pEarthImpl = 0.12
pEarthImpl_unc = 0.04
massBelt = 4e-4
massLV = 4.86e-3
massLV_unc = 1.63e-3
mAsteroids = 607

#nImpl = pImplSD * mAsteroids/massBelt * massLV/pEarthImpl
ImplProb = nImpSim/nSdSim * sd_to_total
ImplProb_unc2 = (1./nSdSim * sd_to_total*nImpSim_unc)**2 + \
                (nImpSim/nSdSim * sd_to_total_unc)**2
print("Implantation probability =", "{:.2E} {:.2E}".format(ImplProb, ImplProb_unc2**0.5))

nImpl = nImpSim/nSdSim * sd_to_total * mAsteroids/massBelt * massLV/pEarthImpl

nImpl_unc2 = (1/nSdSim * sd_to_total * mAsteroids/massBelt * massLV/pEarthImpl * nImpSim_unc)**2 + \
            (nImpSim/nSdSim * 1 * mAsteroids/massBelt * massLV/pEarthImpl * sd_to_total_unc)**2 + \
            (nImpSim/nSdSim * sd_to_total * mAsteroids/massBelt * 1./pEarthImpl * massLV_unc)**2 + \
            (nImpSim/nSdSim * sd_to_total * mAsteroids/massBelt * massLV/pEarthImpl**2 * pEarthImpl_unc)**2

print("Implanted planetesimals =", nImpl, nImpl_unc2**0.5)
